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^•ZJ . We study numerically and analytically a simple ofF-lattice model of scalar harmonic vibrations 

' by means of Euclidean random matrix theory. Since the spectrum of this model shares the most 

puzzling spectral features with the high-frequency domain of glasses (non-Rayleigh broadening of 
the Brillouin peak, boson peak and secondary peak), the Euclidean random matrix theory provide 
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' a single and fairly simple theoretical framework to their explanation. 
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I. INTRODUCTION 



■ The high-frequency (w > IThz) region in vibrational spectra of amorphous systems is related to density fluctuations 
whose size is comparable to the typical distance between particles. Whereas in ordered systems those excitations 
(phonons) persist up to momenta of about the Debye momentum, the fate of excitations of microscopical size in 
^-H ' disordered systems is still a quite puzzling issue, both from the theoretical and the experimental point of view. 

, Recent high-resolution inelastic X-ray scattering (IXS) and neutron scattering techniques have made accessible to 
experiment this region [Q, ||, ^, IJ, |[ ||, |^, ||, |9[ |lO[ |ll|, [ij, 13, |4j where the exchanged external momentum p is 



^2 ! comparable to po, namely the momentum where the static structure factor has its first maximum. A related issue 
is the influence of the high-frequency dynamics of glasses on specific heat experiments in the 1 K-lOO K range. A 
number of facts have emerged from the experiments: 



1. The dynamic structure factor S{p,uj) has a Brillouin-like peak for momenta up to p/po ~ 0.5. This inelastic 
peak is due to the interaction of the external photon (or neutron) with some excitation of the system. A 
very controversial issue is the propagating nature of those excitations . Furthermore a secondary peak 
at frequencies smaller than the Brillouin one develops for larger momenta |l^, |l^, becoming dominant for 
p/po ~ 0.5. 
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' 2. The Brillouin peak width T dependence on the momentum p has been described by means of the following 

scaling law: 

Focp". (1) 

The broadening must generally be ascribed to the interaction among the excitations and the disorder of the 
system. However the details are still unexplained, two remarkable facts indeed arising. First, there seems to be 
quite a general agreement about the fact that F is not affected by changes in the temperature, ruling out the 
hydrodynamic explanation for the broadening. Second, at small enough momenta (e.g. light scattering) a is 
undoubtedly ~ 2, ruling out also Rayleigh scattering. 

3. Beyond the well-known anomalous linear dependence on T at the lowest temperatures of the specific heat (re- 
produced by the two-level model |l^), another deviation from the Debye law at slightly higher temperatures 
is found. The specific heat divided by shows a peak near 10 K ||l^. This peak points to the existence of 
an excess of vibrational states at frequencies ui ~ 1 THz, which shows as a peak in the plot of g{Lu)/uP' {g{uj) 
is the vibrational density of states, as obtained by Raman or inelastic neutron scattering). Since the scattering 
intensity at the peak scales in temperature with Bose statistics, the peak has become known as the Boson peak 
(BP). Two features of the BP should be remarked. First, the peak frequency wbp is several times smaller than 
any natural frequency scale, like the Debye frequency or the band edge. Moreover, where good data for the 
dispersion relation oj{p) (determined from the position of the Brillouin peak of the dynamic structure factor, 
which probes only longitudinal modes) are available, it has been checked that w(p), is still a linear function of 
the momentum p at ujbp ^, ^, |l^. Second, although less experimental data are available, it seems to be the 
rule that uj-bp shifts to lower frequencies on heating [l^, |l^ (except for silica ||ll|, |l^). Yet, in silica the 
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BP evolution upon increasing the density has been studied experimentally [ po| and in simulations pl| | . In close 
agreement with our most recent theoretical results it was found that ujbp shifts to larger frequencies and 
the BP loses intensity when the density grows. The BP has also been identified with the secondary peak in the 
dynamic structure factor at high exchanged momentum, both experimentally [ |T6| and in simulations JTzt . 

A number of basic insights on the spectral properties of glasses have been obtained by means of molecular dynamics 
simulations ||2^, Q on systems such as argon, silica and water. Among others, let us point out 

the fact that, in the glass phase, the high frequency dynamics can be understood in the framework of the harmonic 
approximation (see however the work on soft-potentials pl[ for a dissenting view). 

We still lack an universally agreed upon theoretical interpretation of the propagation of phonons in topologically 
disordered systems. The BP has been variously interpreted as arising from mixing of longitudinal and transverse 
modes ||l^, |l7, 25, 3^, hybridization of optical and acoustic modes [Q, a combination of level repulsion from disorder 
and van Hove singularities j2^ , an associated mechanical instability j2^, ^ , the presence of a loffe-Regel crossover 
at ojbp [§, or from the scattering of sound-waves with localized anharmonic vibrations ||lj, |l|. 

Here we would like to describe a theoretical approach relying on the study of the statistical properties of Euclidean 
Random Matrices (ERM) providing a coherent framework which explains the broadening of the Brillouin peak 
and the origin both of the Boson and of the secondary peak. 

The organization of the paper is as follows. Sec. O de fines our model and gives technical some details on the 
perturbative solution and resummation. Next, sections]! H| and IV discuss the more general theoretical results, while 
in sec. ^ a special interparticle potential is studied, comparing theoretical and simulation results. Sec. VI summarizes 
our conclusions. 



II. REAL SYSTEMS AND ERM THEORY 



We are considering the spectral properties of a system whose particles oscillate around disordered positions. For 
the sake of simplicity we shall consider only the case the scalar case, where all the displacements are collinear. In 
such a system the spring constants are deterministic function, depending only on the distances of the particles rest 
positions, the Hamiltonian being: 

i7=i^(/..M,^,0, (2) 

where M is a matrix of the form: 

N 

Mij = dij^f{xi-Xk)~ f{xi-Xj) , i,j = 1,2 ... ,N . (3) 
fc=i 

If we think of the glass as a disordered harmonic solid, the matrix M is nothing but the Hessian matrix of the system 
where the function /(r) is the second derivative of the pair potential. 

Hence the disorder in the interactions is due to topological reasons, namely the disordered position of particles. As 
in those topologically disordered systems one cannot split the interactions (spring constants) in an ordered part plus 
a disorder-dependent correction, a new theoretical framework involving the so-called Euclidean Random Matrices had 
to be introduced iQ. The quantity to be directly compared with scattering experiments is the dynamic structure 
factor S{p,oj), which at one-excitation level can be computed by means of ||35|| : 

^(1) (p, W) = _ ^^bTp^ ^.^ ^^^^ ^2 ^ (4^ 

u)n 77^0+ 

having introduced the resolvent: 

G(p, z) ^ 1 ^ exp[zp.(=r7-<^)][(z-M)-i],, (5) 

The overline has the meaning of average over the disordered equilibrium positions. As a matter of fact, we are 
assuming that the vibrational spectra of amorphous systems are .self-averaging quantity, hence they should not depend 
on the given configuration of disorder. 
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A straightforward method to compute the resolvent is to write it as the sum of a geometric series, whose i?th 
element is: 



1 e^^^'-oL^.Y^fi' 

ko,ki...kR \ zi 

\ ZR 



Xz.)~f{xkR^,-Xk^)je-^P^''- (6) 

Letting f{p) be the Fourier transform of f{r), without any loss of generality the resolvent can be written: 

G{p,z) = . / 7, (7) 

z - e{p) - z) 

e{p) = p[/(0) -/>)]. (8) 

where the self-energy S describes the interactions between the phonons and the disorder. 

The simplest case arises when all the particles are randomly placed without any correlation. In that situation, using 
P[x\ = {\/V)^ as the probability distribution of quenched variables, it is possible to build a perturbative approach 
where the leading term is exact when the density p in infinite and the next-to leading terms are proportional to powers 
of 1/p p5[ |. Furthermore by means of well-known resummation techniques, the 1/p correction has been exploited in 
order to obtain an integral equation for the self-energy [p7| : 



pj [^T^r 



)D z-e{q)-nq,z) 



(9) 



whose solution takes into account a given class of terms of any order. However, the quenched position in glasses (and 
amorphous systems in general) are highly correlated. Let us sketch a simple approximation dealing with the correlated 
case. Let g'^^\xi^ . . . ,xr) be the i?-points correlation function related to an arbitrary probability distribution P[x\: 

N 

Jl,--- JH 

Hence, the average on the position of the particles is: 

R 



(11) 



Although the computation using the full correlation function would be exceedingly difficult, some progress can be 
made by using the so called superposition approximation: 

g{xi . . . XR+i) = g{xi - X2)g{x2 - x^) . . . g{xR - xr+i) (12) 

where the pair correlation function is used to take into account the correlation of the position of the particles. The 
superposition approximation can be embedded in our calculation if we make the substitution [ p5| : 

/(r)^ .g(r)/(r) (13) 

This is rather important, because, for typical applications, the function /, being badly divergent at short distances, 
does not have a Fourier transform. On the other hand, the function g{r) typically tends to zero at the origin 
exponentially, thus taking care of the algebraic divergence of f{r). 

Let us point out that in the one-excitation approximation the following relationship between the DOS and the 
dynamic structure factor holds ||35[| : 

^'H-r^^^'Hf-^^,^) (14) 

As at very high frequencies the one-excitation approximation does not hold, and many-excitations contributions should 
be taken into account, the reliability of Eq. ( [l4| ) in describing real systems is a very interesting matter. 
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III. BRILLOUIN PEAK 



The following connections between the main features of the dynamical structure factor and the self-energy are 
established: 

• The 'bare' dispersion relation e{p), which would give the position of the peak in the elastic medium limit, is 
renormalized by the real part of the self-energy z). This gives a;''™°''™(p), the position of the maximum 
of the structure factor in the frequency domain. Let us note that uj'^'^^°™[p) is certainly linear for small p, as 
expected. 

• The imaginary part z) computed at the position of the peak uj = gives the width, T{p), of the 
S'-^^p^Lj) by means of: 

= tj'-'="°™(p)r(p) (15) 

Here we want to show that eq. (^ provides a model-independent derivation of the exponent a in the scaling law (|^). 
Indeed, the large q contribution to the imaginary part of the integral in equation ^ is, because of ( [T^ ) 

E^'(p,z) = -^pgA(A) J ^(^f[q)-f[p-q))\ (16) 

where g\{\) is the density of states in the domain of eigenvalues (A — uj'^,gx{uj^) = ^^^)- If the spectrum is Debye- 
like we have g\{X) oc A"'^, and it is straightforward to show that ( p^ ) is proportional to uj^™°^"^ (p) p'^ . Then the 
relation ( [T5| ) implies the scaling: 

ro(p)cxp2 (17) 

irrespective of the function f{r). Clearly this is only the large q contribution to the integral, but it has been shown 
that it indeed controls the peak width at small p . Hence the ERM theory yields the correct asymptotic behavior 
at very low momenta of the broadening of the peak. 



IV. THE BOSON PEAK 



Exploiting the relation (14), it is possible to obtain from (^ an integral equation even for the DOS. As a matter of 



fact, defining G{z) = G{p = oo, z), the DOS turns out to be 

2w 



g{u;)^ ImQiuj-' + iO+). (18) 



TT 



G being the solution of the following integral equation: 



1 £ _ /(O) _ AGiz) - f-fl-p{q)G{q, z), (19) 



pQiz) p J {2tt 

where A — {2tt)^^ Jj{q) d^l- With this equation, one needs to know the resolvent at all q to obtain the DOS, due 
to the last term in the r.h.s. This can be done by solving numerically the self-consistent equation of ref. but 
here we perform an approximate analysis, which is more illuminating. The crudest approximation is to neglect this 
term, in which case Eq. ^is quadratic in Q, and one easily finds a semicircular DOS, with center at uP' = pf{0) and 
radius 2^pA. This spectrum is the glass analogue of a van-Hove singularity: indeed, when p — > oo, the spectrum is 
made of plane waves, with dispersion relation ui'^{p) = p{f{0) — f(p))j P^ , ^ (a continuous elastic medium). uj{p) 
saturates for large p at = pf{0), yielding an enormous pile- up of states which causes the DOS to be concentrated 
at this value ||34| . At finite p, density fluctuations of the x^"^ are present which act as a perturbation that splits 
this degeneracy, and yields the semicircular part of the spectrum at high frequency. But the semicircular spectrum 
misses the Debye part, and a better approximation is needed. So we substitute G in the last term of the r.h.s. by 
the resolvent of the continuum elastic medium Go{z,p) — {z — tj^(p))^^. This is reasonable because the /^(g) factor 
makes low momenta dominate the integral, and due to translational invariance G{z,p) « Go{z,p) in this region 
We shall be looking at small ui, so to a good approximation 

'''^ nq)G,iq,z).-lB-.d^., (20) 
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where the sound velocity is c = y pf "{0)/2 and B > 0. Then Eq. ( |19| ) is again quadratic in and can be solved to 
give 



2pA y ^ [cj2 -p/(0) + S + zp/2(0V/(W)]2y 

We have two limiting cases. At high densities and low frequencies (p/(0) ^ w^, B, 2^/pA), i.e. when the semicircular 
part of the DOS does not reach low frequencies, the square root can be Taylor-expanded, and one gets 

- (22) 

which is precisely Debye's law. At small densities, on the other hand, the center of the semicircle (which is at 
Lu'^ = pf{0) — B) starts to be comparable to its radius (oc y^), meaning that the states in the semicircle hybridize 
with the sound waves. Mathematically, the instability arises when G{0) develops an imaginary part. This can only 
come from the square root in Eq. ([2l|), and it will happen for p < p^, with p^ fixed by the condition 

2^/Ai, + B = pJiO). (23) 

Now when p^ Pc and cj <C ijJ* = 2Trc^y/ pcA/ (pc/^(0)), the square root in Eq. ( pT| ) behaves as 



^D{p- p,)-iuj/u;*, (24) 
with D a positive constant. Here we can distinguish two regimes: 

• When Lu* D{p— pc) ^ w ^ w*, the imaginary part of Q is proportional to ^JZj^ and thus the DOS is g{uj) oc ur"/"^. 

• When uj <C uo*D{p — pc) <C w*, we have an imaginary part (x w, and g{uj) = uP' j [lo* \/ PcAD{p — Pc))- So the 
DOS is Debye-like, but with a very large prefactor, basically unrelated to the speed of sound. 

It is therefore natural to identify wbp with uj*D{p — pc), that can indeed be arbitrarily small upon approaching 
the instability. Notice that the mechanical instability is a kind of phase transition, for which the order parameter is 
— ImC/(0). From Eq. ( |2^ ) we see that this order parameter behaves as {pc — p)^ , with f3 = 1/2, like in mean-field 
theories. 

Since the behavior of the propagator at high momentum does not strongly affect the dispersion relation [ p7| , we do 
not expect deviations from a linear dispersion relation at lobp- This has been checked either numerically or solving 
numerically the self-consistency equation given a particular choice for the function /(r) (see the numerical results 
section below). 



V. AN EXAMPLE: THE NEARLY-GAUSSIAN CASE 



As previously stated, the model is completely defined when / and the distribution of the x are chosen. As discussed 
in sec. (|^, a good approximation taking into account their correlations is the replacement of the spring constants 
with effective ones. That amounts to taking /(r) = g{r)v" {r) instead of v" {r). 

As a model for the effective spring constants, we shall consider the family of functions 

/„(r) = (l-arVa2)e-^^/(2-^), (25) 

where < a < 0.2 (the upper bound has to be imposed in order to guarantee a positive sound velocity). When a = 
(Gaussian case), the Hessian is strictly positive. When a > 0, we have a stable elastic solid at high densities, while 
at low enough densities, typical interparticle distances will be large enough to allow negative eigenvalues (imaginary 
frequencies). Therefore the density p controls the appearance of a mechanical instability in this model. The coun- 
terpart in real glasses of that instability is expected to be the mode-coupling transition. As a matter of fact, that 
dynamical transition marks a change of the local topological properties of the potential energy landscape . As a 
consequence, there is a transition from a region where the short time dynamics is ruled by a positive Hessian (high 
density) to one where an extensive number of negative eigenvalues is found (low density). 

The choice ( p5[ ) may seem an oversimplification, too distant from any realistic case. However, this is not actually so, 
at least for small momenta. For the sake of comparison we can see in figure (^ the Fourier transform of the function 
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g{r)v"{r) for Argon at very low temperature (~ lOK) together with that of ( P5|) in the purely Gaussian case. Since 
the Fourier transform of our force decreases an order of magnitude by pq = 2/a, we shall take this as our po during 
the following discussion, and a will be our unit of length. 

Choosing two different values of a = 0, 0.1, we have both computed the DOS of this model numerically using the 
method of moments with a box of side L — 128 a (more than 5 x 10^ particles), which allows to reconstruct the 
spectrum up to very low frequencies, and numerically solved equation for several values of p, thus obtaining the 
structure factor and the density of states. 



A. Theory vs. Simulations: Dynamic structure factor 

In figure (|^) we show both the S'(p, w) for several values of the momentum (top) and g(cj)/w^ (bottom) as obtained 
from equation (^) at p = 1, together with the results from numerical simulations for the model (|2^) with a = 0. 

Let us note that a very good agreement with the numerical data is achieved. We also found that the agreement is 
still satisfactorily for densities down to p « 0.6. Notice that even for p ~ 1, the cactus resummation fails to reproduce 
the exponential decay at large frequencies of the density of states. 

Finally, let us look at the scaling of the width of the peak in the frequency domain. In the inset of figure (^) we 
plot r(p), obtained by means of ( [l5|) . As expected, the scaling is found for very small momenta, which crosses over 
to a region where a simple law as (bj) is not suitable to describe the real behavior of the system. 

Note that the region where the p scaling is actually found, i.e. p/po < 0.1 is quite different from the region explored 
by X-rays and neutrons scattering experiments, which rather spans the momenta 0.1 < p/pa < 0.5. It is worthwhile 



to note that the same conclusion can be drawn using MCT for hard spheres |26 



B. Theory vs. Simulations: The Boson peak 

In sec. ( |vl ) we have shown that the main result of the ERM theory, namely the integral equation (||), predicts 
the rising of a peak of the function g[ijj)/uj^ in the low frequency region when at low enough densities the system 
approaches a kind of phase transition where negative eigenvalues (imaginary frequencies) begin to appear. 

This is what happens indeed in our model (|25| ) as shown by the DOS obtained with the method of moments. In 
Fig. (^) (top) we show the DOS divided by at several densities for a = 0.1. All these densities are well above the 
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FIG. 2: (Top) Dynamic structure factor for different value of the momentum, obtained both solving the integral equation 
derived from the cactus resummation (theory) and using the method of moments (simulations) in order to get the spectra of a 
vibrational off-lattice model with spring constants given by (|2^) with a = 0. (For lower momenta, the comparison cannot be 
done due to finite volume effects bq]). (Bottom) The comparison between the theoretical and experimental DOS (?(a;) divided 
by the Debye contribution lo'^ for the same model. 
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FIG. 3: The scaling of the broadening a of the Brillouin peak for the a = model as a function of the momentum p. (Inset) 
The large q contribution Fo as obtained by means of eqs. ( |l5| ) ( |l^ ) with respect to the total broadening P. 
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critical density, which for this model is difficult to locate. Anyhow, at p = O.OScr"^ imaginary frequencies are clearly 
found. As predicted, a peak in g{uo)/Lo^ arises on approaching the instability. As density is reduced, the peak grows 
(relative to the Debye value, also plotted) and moves to lower frequencies. The maximum of the peak lies always 
in the linear region of the dispersion relation (Fig. ^, bottom). Since decreasing density plays the role of increasing 
temperature in our model, this peak reproduces the experimental features of the BP. 
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FIG. 4: (Top) Raising of Boson peak in the nearly-Gaussian [a — 0.1) case when the density is lowered down to the critical 
density pc. The BP grows and shifts towards lower frequencies upon decreasing the density. (Bottom) g{uj)/Lj^ and uj{p) for 
p = 0.25/a~^ show that the BP Ues in the region of frequencies where the relation dispersion is still linear. 
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FIG. 5: Raising of Boson peak in the purely Gaussian (a = 0) case. 
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Interestingly enough, the BP even shows up in the purely Gaussian model (see Fig. (|^)), which is stable at all 
densities because a — 0. This case can be thought as the situation where the critical density is pc = 0. 

Since the rising of the BP seems to be a low-density feature of vibrational spectra and the eq. (j^ is best suited 
to describe correctly the high-density region, only a qualitative agreement among the theory and numerics can be 
achieved for the model (^). 

This can be checked looking at Fig. (^), where we show the DOS obtained from the self-consistent G{q, z), and the 
behavior of — Im Q (0) with density. 




0.50 0.51 p 0.52 0.53 



FIG. 6: Top: Reduced DOS divided by for a = 0.1 at several densities, from the numerical solution of the self-consistent 
equation of |^^. Bottom: Order parameter (— ImC/(0)) of the mechanical instability phase transition vs. density (points). Solid 
line is a fit to the predicted behavior A(pc — p)^^^ the fitting parameters being A and pc- 

As stated before, the analytical solution compares reasonably well with the numerical solution (Fig. ^) on a quali- 
tative level, since it greatly overestimates pc- A still open point is the agreement about the critical exponent f3. 

Finally, in the fig. (j?]) there is the numerical confirmation the the relation ( p^ among the DOS g{uj) and infinite 
momentum limit of the dynamic structure factor S'^^-'(p, (at one-excitation level) actually holds. Therefore we 
conclude that the secondary-peak identified in the dynamic structure factor, is nothing but the Boson-Peak. 

VI. CONCLUSIONS 

In summary, we have applied the Euclidean random matrix approach to the study of the high-frequency excitations of 
glassy systems. The main mathematical result is a non-linear integral equation whose solution provides the dynamic 
structure factor (at the one excitation approximation) and the density of states. Our approach has allowed us to 
theoretically confirm, in a model-independent way, that the Brillouin peak of the dynamic structure factor corresponds 
to the propagation of sound-waves in the glass. At small momentum the width of the Brillouin peak scales as p^. At 
higher p (actually, in the experimentally relevant range) a more complicated law is found, but as shown in Fig. |^, the 
local logarithmic derivative is always smaller than that corresponding to Rayleigh's p'^ law. At still larger exchanged- 
momentum, the dynamic structure factor starts to collapse onto the DOS. The frequency and width of the Brillouin 
peak start to be momentum independent. The experimentally identified secondary peak in the dynamic structure 
factor corresponds to the Boson peak in an intermediate p regime in which the Brillouin peak width and position 
stills depends on exchanged momentum. 
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FIG. 7: The function g{ijj)/uj^ (labelled as DOS) and the function S{p,ui) / {KbTp^) at different values of the momentum p 
given hy p — 2-KN/a for the q = 0.1 model a,t p — 0.25/ a^^ . 

Our theoretical approach is also able to cope with the Boson peak. In the family of models that we consider, there 
is a mechanical instability transition controlled by the density, as signalled by the presence of imaginary frequencies. 
The vibrational density of states of our models contains a BP which is the precursor of the instability transition. 
The BP in our model shares the main features of the experimental BP: it appears for frequencies in the linear part 
of the dispersion relation and it shifts towards arbitrarily low frequencies on approaching a mechanical instability. 
We also reproduce qualitatively the behavior of the silica BP when the density changes ^ (the detailed theory 
of the temperature evolution of the silica BP should consider its negative thermal dilatation coefficient). The BP 
is built from the hybridization of sound waves with high frequency modes (extended but non propagating) that get 
softer upon approaching the instability. The analogous of our instability transition in nature is the topological phase 
transition |36| that underlies the dynamic crossover at the Mode Coupling temperature of real glasses The 
precise nature of the high-frequency modes that hybridize with the sound-waves is most likely material dependent 
and non-universal: they could be transverse ^ ^ ^ , optical , or even longitudinal modes as in our model. 
We believe however that the basic mechanism for the formation of the BP uncovered in our model is common to most 
(if not all) structural glasses. Yet all real glasses do have transverse excitations, and one could ask about generic new 
features introduced by these modes. Work is currently in progress to extend the ERM approach to include transverse 
displacements. 
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